# =============================================================== #
# Reproducibility of:                                             #
# Drivers and Trends for the Equality of Opportunity for Sexual   #
# and Gender Minorities: A Panel Approach Equality of Opportunity #
# for Sexual and Gender Minorities 2024                           #
#                                                                 #
# Code written by Omar Alburqueque and reviewed by Paola Ballon   #
# Contact: oalburquequechav@worldbank.org, pballon@worldbank.org  #
# =============================================================== #

# ------------- #
# Tables A4-A5  #
# ------------- #

library(writexl)

# Locate project root.
# This works whether the script is run from the project root or from /scripts.
wd <- normalizePath(getwd(), winslash = "/", mustWork = TRUE)

if (dir.exists(file.path(wd, "outputs")) &&
    dir.exists(file.path(wd, "intermediate files"))) {
  project_root <- wd
} else if (dir.exists(file.path(dirname(wd), "outputs")) &&
           dir.exists(file.path(dirname(wd), "intermediate files"))) {
  project_root <- dirname(wd)
} else {
  stop("Project root not found. Run this script from the project root or from the scripts folder.")
}

intermediate_dir <- file.path(project_root, "intermediate files")
outputs_dir      <- file.path(project_root, "outputs")

# --------------------------------------------------------------- #
# Function to parse first-stage output from Stata xtivreg logs     #
# --------------------------------------------------------------- #

parse_xtivreg_firststage <- function(path_log, colname) {
  
  if (!file.exists(path_log)) {
    stop(paste("Log file not found:", path_log))
  }
  
  txt <- readLines(path_log, warn = FALSE)
  
  # Locate first-stage section
  i_first <- grep("^First-stage G2SLS regression", txt)[1]
  if (is.na(i_first)) {
    stop(paste("First-stage section not found in:", path_log))
  }
  
  # Number of observations
  i_no <- grep("Number of obs", txt)
  i_no <- i_no[i_no > i_first][1]
  
  n_obs <- if (!is.na(i_no)) {
    as.integer(gsub("[^0-9]", "", txt[i_no]))
  } else {
    NA_integer_
  }
  
  # Number of countries/clusters
  i_clust <- grep("adjusted\\s+for\\s+[0-9,]+\\s+clusters",
                  txt, ignore.case = TRUE)
  i_clust <- i_clust[i_clust > i_first][1]
  
  n_ctry <- if (!is.na(i_clust)) {
    as.integer(gsub(",", "", sub(".*for\\s+([0-9,]+)\\s+clusters.*",
                                 "\\1", txt[i_clust])))
  } else {
    NA_integer_
  }
  
  # Locate coefficient table
  i_header <- grep("\\|.*(Coef\\.|Coefficient)", txt, ignore.case = TRUE)
  i_header <- i_header[i_header > i_first][1]
  
  if (is.na(i_header)) {
    stop(paste("Coefficient header not found in:", path_log))
  }
  
  i_plus <- which(grepl("^-+\\+[-]+", txt))
  i_plus <- i_plus[i_plus > i_header][1]
  
  if (is.na(i_plus)) {
    stop(paste("Start of coefficient table not found in:", path_log))
  }
  
  i_plain <- which(grepl("^-{5,}\\s*$", txt))
  i_plain <- i_plain[i_plain > i_plus][1]
  
  if (is.na(i_plain)) {
    stop(paste("End of coefficient table not found in:", path_log))
  }
  
  lines <- txt[(i_plus + 1):(i_plain - 1)]
  lines <- lines[nzchar(trimws(lines))]
  lines <- lines[grepl("\\|", lines)]
  
  lines_clean <- gsub("\\s+", " ", trimws(lines))
  parts <- strsplit(lines_clean, "\\s*\\|\\s*")
  
  left  <- vapply(parts, `[`, character(1), 1)
  right <- vapply(parts, `[`, character(1), 2)
  
  nums_list <- strsplit(right, " ")
  
  nums_mat <- t(sapply(nums_list, function(x) {
    x <- x[x != ""]
    x <- x[1:6]
    as.numeric(x)
  }))
  
  colnames(nums_mat) <- c("coef", "se", "z", "p", "ci_lo", "ci_hi")
  
  df <- data.frame(
    variable = left,
    nums_mat,
    stringsAsFactors = FALSE,
    row.names = NULL
  )
  
  # Keep and order the variables needed in Tables A4 and A5
  expected_vars <- c(
    "ipgi",
    "pogi",
    "vdem_l20",
    "ltrad_1",
    "ltrad_2",
    "figidj",
    "trgidj",
    "rnna",
    "_cons"
  )
  
  missing_vars <- setdiff(expected_vars, df$variable)
  
  if (length(missing_vars) > 0) {
    stop(paste(
      "Missing variables in first-stage log:",
      paste(missing_vars, collapse = ", "),
      "in",
      path_log
    ))
  }
  
  df <- df[match(expected_vars, df$variable), ]
  
  variable_labels <- c(
    ipgi     = "Interpersonal Globalization Index",
    pogi     = "Political Globalization Index",
    vdem_l20 = "Participatory Democracy Index (lag)",
    ltrad_1  = "Colonized, common-law system",
    ltrad_2  = "Colonized, civil-law system",
    figidj   = "Financial Globalization Index (de jure)",
    trgidj   = "Trade Globalization Index (de jure)",
    rnna     = "Log of capital stock at constant 2017 national prices",
    `_cons`  = "Constant"
  )
  
  stars <- ifelse(df$p < 0.01, "***",
                  ifelse(df$p < 0.05, "**",
                         ifelse(df$p < 0.10, "*", "")))
  
  coef_str <- paste0(sprintf("%.2f", df$coef), stars)
  se_str   <- paste0("(", sprintf("%.2f", df$se), ")")
  
  value_col <- as.vector(rbind(coef_str, se_str))
  
  var_col <- rep(variable_labels[df$variable], each = 2)
  var_col[seq(2, length(var_col), by = 2)] <- ""
  
  value_col <- c(
    value_col,
    if (!is.na(n_obs)) format(n_obs, big.mark = ",", scientific = FALSE) else "",
    if (!is.na(n_ctry)) as.character(n_ctry) else ""
  )
  
  var_col <- c(var_col, "Obs.", "Countries")
  
  out <- data.frame(
    Variable = var_col,
    value_col,
    stringsAsFactors = FALSE
  )
  
  names(out)[2] <- colname
  
  out
}

# -------- #
# Table A4 #
# -------- #

tableA4 <- parse_xtivreg_firststage(
  file.path(intermediate_dir, "xtivreg_firststage.log"),
  "(4) RE IV + Colonial condition-legal system"
)

write_xlsx(tableA4, file.path(outputs_dir, "TableA4.xlsx"))

# -------- #
# Table A5 #
# -------- #

dim1 <- parse_xtivreg_firststage(
  file.path(intermediate_dir, "xtivreg_firststage_dim1.log"),
  "Decriminalization"
)

dim2 <- parse_xtivreg_firststage(
  file.path(intermediate_dir, "xtivreg_firststage_dim2.log"),
  "Access to Education"
)

dim3 <- parse_xtivreg_firststage(
  file.path(intermediate_dir, "xtivreg_firststage_dim3.log"),
  "Access to Labor Markets"
)

dim4 <- parse_xtivreg_firststage(
  file.path(intermediate_dir, "xtivreg_firststage_dim4.log"),
  "Access to Services and Social Protection"
)

dim5 <- parse_xtivreg_firststage(
  file.path(intermediate_dir, "xtivreg_firststage_dim5.log"),
  "Civil and Political Inclusion"
)

dim6 <- parse_xtivreg_firststage(
  file.path(intermediate_dir, "xtivreg_firststage_dim6.log"),
  "Protection from Hate Crimes"
)

tableA5 <- data.frame(
  Variable = dim1$Variable,
  Decriminalization = dim1[[2]],
  `Access to Education` = dim2[[2]],
  `Access to Labor Markets` = dim3[[2]],
  `Access to Services and Social Protection` = dim4[[2]],
  `Civil and Political Inclusion` = dim5[[2]],
  `Protection from Hate Crimes` = dim6[[2]],
  check.names = FALSE,
  stringsAsFactors = FALSE
)

# Add Sargan-Hansen rows from Stata diagnostics
overid_path <- file.path(intermediate_dir, "iv_overid_diagnostics.csv")

if (!file.exists(overid_path)) {
  stop(paste("Diagnostics file not found:", overid_path))
}

overid <- read.csv(overid_path, stringsAsFactors = FALSE)

get_overid_value <- function(outcome_name, statistic_name) {
  
  row <- overid[
    overid$table == "Table A5" &
      overid$outcome == outcome_name,
  ]
  
  if (nrow(row) == 0) {
    return("")
  }
  
  val <- row[[statistic_name]][1]
  
  if (is.na(val)) {
    return("")
  }
  
  if (statistic_name == "sargan_hansen_chi2") {
    return(sprintf("%.3f", val))
  }
  
  if (statistic_name == "sargan_hansen_p") {
    return(sprintf("%.4f", val))
  }
  
  return(as.character(val))
}

sargan_chi2_row <- data.frame(
  Variable = "Sargan-Hansen \u03c7\u00b2(2)",
  Decriminalization = get_overid_value("Decriminalization", "sargan_hansen_chi2"),
  `Access to Education` = get_overid_value("Access to Education", "sargan_hansen_chi2"),
  `Access to Labor Markets` = get_overid_value("Access to Labor Markets", "sargan_hansen_chi2"),
  `Access to Services and Social Protection` = get_overid_value("Access to Services and Social Protection", "sargan_hansen_chi2"),
  `Civil and Political Inclusion` = get_overid_value("Civil and Political Inclusion", "sargan_hansen_chi2"),
  `Protection from Hate Crimes` = get_overid_value("Protection from Hate Crimes", "sargan_hansen_chi2"),
  check.names = FALSE,
  stringsAsFactors = FALSE
)

sargan_p_row <- data.frame(
  Variable = "Sargan-Hansen p-value",
  Decriminalization = get_overid_value("Decriminalization", "sargan_hansen_p"),
  `Access to Education` = get_overid_value("Access to Education", "sargan_hansen_p"),
  `Access to Labor Markets` = get_overid_value("Access to Labor Markets", "sargan_hansen_p"),
  `Access to Services and Social Protection` = get_overid_value("Access to Services and Social Protection", "sargan_hansen_p"),
  `Civil and Political Inclusion` = get_overid_value("Civil and Political Inclusion", "sargan_hansen_p"),
  `Protection from Hate Crimes` = get_overid_value("Protection from Hate Crimes", "sargan_hansen_p"),
  check.names = FALSE,
  stringsAsFactors = FALSE
)

tableA5 <- rbind(tableA5, sargan_chi2_row, sargan_p_row)

write_xlsx(tableA5, file.path(outputs_dir, "TableA5.xlsx"))